clear all;
clc;
%using a Matlab interpolation to get intermediate values for population
%Piecewise Cubic Hermite Interpolating Polynomial (PCHIP)


SN=50; %state numbers (here we do work without DC)
YEARS=42; %making number of years a variable
load ('input1960-2001all.mat');

%reversing the stacking to use code
for o=1:YEARS
x(:,o)=pop(1+SN*(o-1):SN*o);
end;
pop=x;
clear x

for o=1:YEARS
x(:,o)=income(1+SN*(o-1):SN*o);
end;
income=x;
clear x

for o=1:YEARS
x(:,o)=incomepc(1+SN*(o-1):SN*o);
end;
incomepc=x;
clear x

[pop80] = xlsread('pop1980.xls');     % population in millions 
[pop90] = xlsread('pop1990.xls');     % population in millions 
[pop00] = xlsread('pop2000.xls');     % population in millions 

Income=income;
Incomepc=incomepc;

for i=1:SN
x = [1980,1990,2000]; 
y(i,:) = [pop80(i),pop90(i),pop00(i)]; 
t = [1980:2001];
p(i,:) = pchip(x,y(i,:),t);
end;

popforecast80=[pop(:,1:20) p];


%growth rates for Population in sample und out of sample 
            
   for i=2:YEARS
      popgth(:,i-1)=(popforecast80(:,i)-popforecast80(:,i-1))./popforecast80(:,i);
  end 
       clear i;
 popgth80=popgth;        
       
%growth rates of income in sample
 for i=2:21
      incomegth(:,i-1)=(income(:,i)-income(:,(i-1)))./income(:,i);
  end 
       clear i;
            
% Create bivariate normal distribution for population and GDP growth rate          

for i=1:SN   
    pm=mean(popgth(i,1:20));
    ps = std(popgth(i,1:20));
   
    g=mean(incomegth(i,1:20));
    gs=std(incomegth(i,1:20));

    rho= corrcoef([popgth(i,1:20)', incomegth(i,1:20)']);
    
    rho=rho(1,2);
    
    beta=rho*ps*gs/(ps*ps);
            
            for v=1:21;
                incomegth(i,(20+v))=(g-beta*pm+beta*popgth(i,(20+v))) + (gs*sqrt(1-(rho*rho)))*randn(1);
            end;
            
end
clear i v
   % Calculating income Projection for all states
            for i=1:21;
                income(1:SN,(21+i))=income(1:SN,(20+i)).*(ones(SN,1)+incomegth(1:SN,(20+i)));
            end;
            clear i;
            
     
    for i=1:21           
  incomepc(1:SN,(21+i))=income(1:SN,(21+i))./(popforecast80(1:SN,(21+i)));
    end
clear i

for i=1:41           
  incomepcgth(1:SN,i)=(incomepc(1:SN,(i+1))-incomepc(1:SN,(i)))./incomepc(1:SN,(i));
end
incomepcgth80=incomepcgth;
clear i
t=1960:2001;
for i=[2,7,10,24,29,35,41]
 figure(i)
 plot(t,incomepc(i,:),'yellow')
 hold
 plot(t,Incomepc(i,:))
 
end

incomeforecast80=income;
incomepcforecast80=incomepc;

popdensforecast80=(popforecast80(1:SN,:)*1000000)./kron(ones(1,YEARS),area50(1:SN));

yangschneidergth80factor=(ones(SN,21)+incomepcgth80(1:SN,21:41)).*(ones(SN,21)+popgth(1:SN,21:41)).*kron(ones(SN,21),(1-0.01697)*(1+0.00136));